'''
Created on Oct 8, 2010

@author: oabalbin
'''
import sys
import numpy as np
from optparse import OptionParser
from datetime import datetime
from collections import defaultdict, deque

import RNAseq.array as seqarray
import RNAseq.sampleinfo.sample_info as si
import RNAseq.arraytools.subset_genes as sg
import RNAseq.arraytools.copa as ac
import RNAseq.io.translate_gene_names as tn

if __name__ == '__main__':
    
    optionparser = OptionParser("usage: %prog [options] ")
    optionparser.add_option("-f", "--annotFile", dest="annotFile",
                            help="annotation file for all files to use")
    optionparser.add_option("-d", "--outfolder", dest="outfolder",
                            help="output folder")
    optionparser.add_option("-s", "--samplesFile", dest="samplesFile",
                            help="annotation file with the sample information")
    optionparser.add_option("-g", "--genesFolder", dest="genesFolder", default="use_all_genes",
                            help="annotation file with genes of interest")    
    optionparser.add_option("-t", "--tissueName", dest="tissueName", default="all_tissues",
                            help="name of the tissue analyzed")
    optionparser.add_option("-p", "--sampleType", dest="sampleType",
                            help="sample type. could be Cell Line or Tissue")
    

    (options, args) = optionparser.parse_args()
    
    
    t = datetime.now()
    tstamp = t.strftime("%Y_%m_%d_%H_%M")
            
    thisfile = open(options.annotFile)
    #outfile_outliers_list = open(time_stamp_outfiles[0],'w')
    #outfile_outliers_matrix = open(time_stamp_outfiles[1],'w')
    array_dump_file = options.annotFile + '_dump_file'

    ####################
    ##### Parameters
    percentils_list =  [80,85,90,95,98] #[75,80,85,90] #[80,85,90,95,98] 
    mmean=True               # median, mean, sum

    nout2report=10
    expression_cutoff=0.0       # minimum median RPKM of expression
    base_expression_cutoff= 50  # In percentil. It is used to report the expression value of a particular gen in a sample. 
                                # The median usually can be used, however because in many cases median=1, It is maybe desireable to use for e.g the75%. 
    median_shift=1.0
    sdv_factor=0.25             #0.25
    not_report_isoforms=True    # Report gene outlier isoforms
    nrandom_permutations = 10000 # Number of random permutations
    
    
    # create the meme array
    filename=options.annotFile
    matrixfilename=options.annotFile+'memmap'
    header_cols=2
    matobj = seqarray.read_matrix(filename, matrixfilename, header_cols)
    all_samples, all_genes, all_gene_intervals = \
    seqarray.create_matobj_dicts(matobj.cols,matobj.sample_names)
    #print set(matobj.sample_params['diagnosis'])
    cohorts = set(matobj.sample_params['cohort'])
    ####################
    # samples and gene list to use

    if options.genesFolder!='use_all_genes':
        # Get all the gene files in the folder
        gene_list_file = options.genesFolder
        use_all_genes=False
    else:
        use_all_genes=True
        
    ########### Program          
    # subset of genes and samples to consider 
    # Get the subset of genes
    if not use_all_genes:
        thisgenelist = sg.list_of_names(open(gene_list_file))   
        gene_subset, gene_rows, gene_subset_intervals = sg.get_gene_subset_mod(thisgenelist, all_genes, all_gene_intervals)
    else:
        gene_subset, gene_rows, gene_subset_intervals = sg.get_gene_subset_mod(all_genes.keys(), all_genes, all_gene_intervals)
    
    if options.tissueName != 'all_tissues':
        outfiles = ['outfile_outliers_list','outfile_outliers_matrix']
        time_stamp_outfiles=[]
        for name in outfiles: 
            time_stamp_outfiles.append( options.outfolder+tstamp+'_'+options.sampleType+'_'+options.tissueName+'_'+name+'.txt')

        outfile_outliers_list = open(time_stamp_outfiles[0],'w')
        outfile_outliers_matrix = open(time_stamp_outfiles[1],'w')
        
        select_samples_by = [options.sampleType.replace('_',' '),options.tissueName,'PASS']
        sampleInfoFile = open(options.samplesFile)
        sample_info_dict = si.get_sample_info(sampleInfoFile,select_samples_by)
        sample_dict = si.generate_sample_lists(sample_info_dict)
        sample_subset, sample_cols, benign_cols = si.get_sample_columns(sample_dict, all_samples)
        tumor_cols = list(set(sample_subset.keys()).difference(benign_cols))
        cohort_sizes = [len(sample_cols), len(tumor_cols), len(benign_cols)]
        print sample_dict
            
        # create a sub array  using the genelist and the samplelist
        expression_matrix = matobj.exprs[gene_rows, :]
        expression_matrix = expression_matrix[:,sample_cols]
        
        ### Start the copa analysis
        # Returns: outlier list to report, all gene-copa scores, dictionary with gene-outlier sample, matrix of outliers by samples
        outliers_list, gene_copa_score_array, gene_outliersample_dict, matrix_outliers_bysample = \
         ac.copa_analysis(expression_matrix, median_shift, benign_cols, tumor_cols, percentils_list, sdv_factor, mmean)
        
        # Calculate the pvalue for the copa scores
        print "Calculating permutations"
        copa_pvalues = ac.calculate_score_pvalue(gene_copa_score_array, nrandom_permutations)
        
        #### Print
        inputfile_genes = '/data/ucsc_tables/knownGene.txt'
        inputfile_Hugo = '/data/ucsc_tables/kgXref.txt'
        knownGenes = tn.read_known_genes(inputfile_genes)
        knownHugoNames = tn.read_hugo2ucsc_names(inputfile_Hugo)
        
        ac.print_list_of_outliers(outliers_list, outfile_outliers_list, copa_pvalues, gene_subset, gene_subset_intervals, knownHugoNames, cohort_sizes)
        ac.print_matrix_of_outliers(sample_subset, gene_subset, outliers_list, 
                                    matrix_outliers_bysample, outfile_outliers_matrix)
    
    else:
        for tissue in cohorts:
                        
            select_samples_by = [options.sampleType.replace('_',' '),tissue,'PASS']
            sampleInfoFile = open(options.samplesFile)
            sample_info_dict = si.get_sample_info(sampleInfoFile,select_samples_by)
            sample_dict = si.generate_sample_lists(sample_info_dict)
            sample_subset, sample_cols, benign_cols = si.get_sample_columns(sample_dict, all_samples)
            tumor_cols = list(set(sample_subset.keys()).difference(benign_cols))
            cohort_sizes = [len(sample_cols), len(tumor_cols), len(benign_cols)]
            
            print sample_dict
            print "Tissue Name ", tissue, len(sum(sample_dict.values(),[]))
            
            if tissue == 'None' or len(sample_cols) < 5: # len(sum(sample_dict.values(),[])) < 5:
                print "Tissue not analyzed", tissue, len(sample_cols) #len(sum(sample_dict.values(),[]))
                continue
            #continue 
            outfiles = ['outfile_outliers_list','outfile_outliers_matrix']
            time_stamp_outfiles=[]
            for name in outfiles: 
                time_stamp_outfiles.append( options.outfolder+tstamp+'_'+options.sampleType+'_'+tissue+'_'+name+'.txt')
    
            outfile_outliers_list = open(time_stamp_outfiles[0],'w')
            outfile_outliers_matrix = open(time_stamp_outfiles[1],'w')

            
            # create a sub array  using the genelist and the samplelist
            expression_matrix = matobj.exprs[gene_rows, :]
            expression_matrix = expression_matrix[:,sample_cols]
            
            ### Start the copa analysis
            # Returns: outlier list to report, all gene-copa scores, dictionary with gene-outlier sample, matrix of outliers by samples
            outliers_list, gene_copa_score_array, gene_outliersample_dict, matrix_outliers_bysample = \
             ac.copa_analysis(expression_matrix, median_shift, benign_cols, tumor_cols, percentils_list, sdv_factor, mmean)
            
            # Calculate the pvalue for the copa scores
            print "Calculating permutations"
            copa_pvalues = ac.calculate_score_pvalue(gene_copa_score_array, nrandom_permutations)
            
            #### Print
            inputfile_genes = '/data/ucsc_tables/knownGene.txt'
            inputfile_Hugo = '/data/ucsc_tables/kgXref.txt'
            knownGenes = tn.read_known_genes(inputfile_genes)
            knownHugoNames = tn.read_hugo2ucsc_names(inputfile_Hugo)
            
            ac.print_list_of_outliers(outliers_list, outfile_outliers_list, copa_pvalues, gene_subset, gene_subset_intervals, knownHugoNames,\
                                      cohort_sizes)
            ac.print_matrix_of_outliers(sample_subset, gene_subset, outliers_list, 
                                        matrix_outliers_bysample, outfile_outliers_matrix)

            
    
    

    sys.exit(0)
